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ABSTRACT 



We develop an analytical equation of state (EOS) for magnetized fully ionized plasmas, which cover a wide range of temperatures and 
densities, from low-density classical plasmas to relativistic, quantum plasma conditions. This EOS directly applies to calculations of 
structure and evolution of strongly magnetized white dwarfs and neutron stars. We review available analytical and numerical results 
for thermodynamic functions of the nonmagnetized and magnetized Coulomb gases, liquids, and solids. We propose a new analytical 
expression for the free energy of solid Coulomb mixtures. Based on recent numerical results, we construct analytical approximations 
for the thermodynamic functions of harmonic Coulomb crystals in quantizing magnetic fields. The analytical description ensures a 
consistent evaluation of all astrophysically important thermodynamic functions based on the first, second, and mixed derivatives of 
the free energy. Our numerical code for calculation of thermodynamic functions based on these approximations is made publicly 
available. Using this code, we calculate and discuss the effects of electron screening and magnetic quantization on the position of the 
melting point in a range of densities and magnetic fields relevant to white dwarfs and outer envelopes of neutron stars. We consider 
also the thermal and mechanical structure of a magnetar envelope and argue that it can have a frozen surface which covers the liquid 
ocean above the solid crust. 
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Coulomb plasmas, i.e., the fully ionized plasmas whose thermo- 
dynamics is strongly affected by electrostatic interactions, are 
encountered in many physical and astrophysical situations (e.g. 
lFortovl l2009'). Full ionization is reached either at high tempera- 
tures T and low densities p (thermal ionization) or at high densi- 
ties p (pressure ionization). The latter case is typical of the inte- 
rior conditions of low-mass stars, brown dwarfs or giant planets 
jChabrier & Baraffe 2000) as well as the interior and envelope 
conditions of white dwarfs and neutron stars. Coulomb interac- 
tions are crucial for the equation of state (EOS) under such con- 
ditions. In the interior or the envelope of compact objects such 
as white dwarfs and neutron stars, the electrons can be weakly 
or strongly degenerate, the plasma can be in the liquid or solid 
state, the electrons can have various degrees of degeneracy and 
relativism, and the quantum effects on ion motion can be sub- 
stantial. Therefore, a wide-range EOS is needed for calculations 
of the structure and evolution of such stars. 

In a previous work (Potekhin & Chabriel^l2000l |2010[ here- 
after Paper I and Paper II, respectively) we proposed a set of 
analytical expressions for the calculations of the EOSs of the 
Coulomb plasmas without magnetic fields and presented a code 
for thermodynamic functions based of the first, second, and 
mixed derivatives of the analytical Helmholtz free energy F with 
respect to density p and temperature T. This code has been em- 
ployed in astrophysical modeling and adapted for the use in the 
MESA project (Paxton et al. 2011i) . 
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The Bohr - van Leeuwen theorem states that an EOS of 
charged po intlike classical par ticles is not affected by mag- 
netic field dvan LeeuwenI 1 1 92lb . However, magnetic field can 
affect thermodynamic functions through intrinsic magnetic mo- 
ments of particles and by quantization of charged particles mo- 
tion in Landau orbitals (lLandatJl930tlLandau & Lifshitzll 19771) . 
These effects can be important, for example, in magnetic white 
dwarfs whose magnetic fields B c an reach 10^ - 10^ G (e.g., 
IWickramasin^he & Ferrariol |2000[ and references therein) and 
neutr on stars with typical B ~ 10^ - 10'"* G (e.g., Haen sel et aP 
l2007l and references therein). 

In this paper, we systematically consider analytical expres- 
sions for thermodynamic functions of magnetized Coulomb 
plasmas, discuss their validity range, and introduce some prac- 
tical modifications. We take account of analytical and numeri- 
cal results, currently available for various contributions to the 
Helmholtz free energy in quantizing magnetic fields. Taking ad- 
vantage of recently published numerical results (Baiko 2009^), 
we construct an analytical description of the thermodynamic 
functions of harmonic Coulomb crystals in quantizing magnetic 
fields. 

In Sect. |2] we give definitions and simple estimates for 
the plasma parameters that determine different thermodynamic 
regimes. In Sect. |3] we outline the EOS of a nonmagnetized 
Coulomb plasma as the reference case. In Sect. |4]we consider 
the Boltzmann and Fermi gases in quantizing magnetic fields, 
present a general analytical description of their EOS, and sim- 
plify them for several limiting cases. In Sect. |5]we review non- 
ideal contributions to the EOS of a Coulomb liquid in a strong 
magnetic field. In Sect. |6] we derive an analytical approxima- 
tion for the EOS of a strongly magnetized Coulomb crystal. In 
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Sect.|2]we present and discuss examples of thermodynamic func- 
tions for conditions typical of white dwarfs and neutron-star en- 
velopes. The summary is given in Sect. |8] In Appendices we 
give the explicit expressions for thermodynamic functions that 
are used in Sect. |3] 



plasma frequency. Since Aija-i — (Tp/T) y2/T, the importance of 
the quantum effects in strongly coupled plasmas (i.e., at F » 1) 
is determined by parameter 

77 = rp/r = rV37^* 7.832 (z/A)^/r6. (6) 



2. Basic definitions 

2.1. General parameters 

Let He and ni be the electron and ion number densities, A and 
Z the ion mass and charge numbers, respectively. In this paper 
we consider the neutral plasmas (therefore - Z«i) that con- 
tain a single type of ions and include neither positrons (they can 
be described using the same analytical fun ctions as for the elec- 
trons; see, e.g., iBlinnikov et al. 1996; Ti mmes & A mett 1999]), 
nor free neutrons (see Haensel et al.ll2007l for a review). 

The state of a free electron gas is determined by the elec- 
tron number density and temperature T. Instead of «e it is 
convenient to introduce the dimensionless density parameter 
Ts - fle/flo, where oq is the Bohr radius and - (^nriey^^^. 
The parameter can be quickly evaluated from the relations 



1.1723n 



-1/3 
24 



(po/p)'^"', where n24 = ne/10^'*cm ^ and 
po - 2.6752 (A/Z) g cm"-'. The analogous density parameter for 
the ions is Rs = aimiiZef jti^ = 1822.89AZ'^/Vs, where nii is 
the ion mass and ai = {^mii}^^!^ is the ion sphere radius. 

At stellar densities it is convenient to use, instead of r^, the 
nonmagnetic relativity parameter 



X, = ^-p/mec = 1.00884 (p6Z/A)'/^= 0.014005 r\ 



(1) 



where pf - h^^n^n^y^^ is the electron Fermi momentum in 
the absence of a magnetic field, and pe = g cm""*. The 

Fermi energy (without the rest energy) for the electron gas is 
ep = c -y'(mec)^ + (pp)^ - m^c^, and the Fermi temperature 
Tp = e p/kB = T, (jr - 1), where T, = m^c^/kji = 5.93 x 10'' K, 
7r = Vl + -''^r^ ^iid is the Boltzmann constant. A useful mea- 
sure of electron degeneracy is = T/Tp. In the nonrelativistic 
hmit (jcr <K 1), Tp * 1.163 X 10'' r;^ K, and 



= 0.543 rJFe 
where 



22.747 / Z\'/^ 



aeksT 



mil zv 



(2) 



(3) 



In the opposite ultrarelativistic limit (x, » 1), ~ (263 Fg)"'. 

The strength of the Coulomb interaction of nonrelativistic 
ions is characterized by the Coulomb coupling parameter 



r . = FeZ^/\ 

QikYiT 



(4) 



where Te = K. 

Thermal de Broglie wavelengths of free ions and electrons 
are usually defined as 



[niiksT 



2 \l/2 



2nti 

MeksT 



2 \l/2 



(5) 



(though in some publications such definitions differ by a numer- 
ical factor). The quantum effects on ion motion are important 
either at A; > ai or at T <K Tp, where Tp = Hojp/kB is the 

/ , - -.1/2 

ion plasma temperature, and cOp — [Ane riiZ jmA is the ion 



2.2. Magnetic-field parameters 

In the nonrelativistic theory (iLandau & Lifshitzlll977l) . the en- 
ergy of an electron in magnetic field B equals nhu^ + me/??/2, 
where p, is the momentum component along B, - eB/nigC 
is the electron cyclotron frequency, « = «l + 5 - charac- 
terizes a Landau level, cr = ±1 determines the spin projection 
on the field, and hl is the non-negative integer Landau number 
related to the quantization of the kinetic motion transverse to 
the field. In the relativi stic theory (iJohnson & Lippmannlll949l : 
Berestets kii et al. I IT981 . the kinetic energy e of an electron at 
the Landau level n and its longitudinal momentum p, are inter- 
related as 

(99 9 \ 9 

m^c + ItuD^m^n + p A - m^c , (7) 

\Pz\ - Pn{e) - [('«eC + e/cf - {m^cf- - 2mJiUi,n^ ^ . (8) 

The levels n > 1 are double-degenerate with respect to cr. Their 
splitting due to the anomalous magnetic moment of the electron 
is K (mf.c^af /2jT) b at « 1 and ~ {meC^af/2n) \lnb - 1.584]^ 
at » 1 (see ISchwingeJil 9881: ISuh & Mathewsll200lh . which is 
always much smaller than ha>c and is negligible in the compact 
stars. 

Convenient dimensionless parameters that characterize the 
magnetic field in a plasma are the ratios of the electron cyclotron 
energy fiwc to the Hartree unit of energy, to the electron rest 
energy, and to A;Br: 

= h^B/mlce^ = B/Bo, (9) 

where Bn = 2.3505 x lO"* G, 



B 



(10) 



(11) 



mec- * 4.414 X 1013 G" 

where Uf - e^/hc is the fine-structure constant, and 

( = hojJkBT ^ 134.345,2/^6, 

where Bn = B/IO'^ G. The magnetic length = (tic/eBy'^ = 
flo/ ^fym gives a characteristic transverse scale of the electron 
wave function. 

For the ions, the cyclotron energy is hoj^ = Z (111^/1111) hojc, 
and the parameter analogous to ( is 

^, = naj,,/kBT ^ 0.0737 {Z/A)Bn/T6. (12) 

Another important parameter is the ratio of the ion cyclotron 
frequency to the plasma frequency. 



yS = Wei/Wp = * 0.0094 B12/ 



2.3. Free energy and thermodynamic functions 



(13) 



The Helmholtz free energy F of a plasma can be conveniently 
written as 



F = + + Fee + Fa + Fie, 



(14) 



where F^' and F^* denote the ideal free energy of the ions and 
the electrons, and the last three terms represent an excess free en- 
ergy arising from the electron-electron, ion-ion, and ion-electron 
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interactions, respectively. In the nonideal plasmas, correlations 
between any plasma particles depend on all interactions, there- 
fore the separation in Eq. (fT4t is just a question of convenience. 

An important reference case is the model of one-component 
plasma (OCP). In this model, the electrons are replaced by a 
"rigid" (nonpolarizable) background of the uniform charge dis- 
tribution. It is convenient to define Fa as the difference between 
F and f'^' in the OCP model. Still stronger simplification is the 
ion-sphere model, in which the interaction energy in the OCP 
is evaluated as the electrostatic energy of a positive ion in the 
negatively charged sphere of radius oi (Salpeter 1961). The elec- 
tron exchange-correlation term is defined as F^e - F - Ff^ in 
the model of an electron gas without consideration of the ions, 
which are replaced by an uniform positive background to ensure 
the global charge neutrality. Then the ion-electron (electron po- 
larization) contribution Fie is the difference between F and the 
other terms, when interactions between all types of particles are 
taken into account. 

The pressure P, the internal energy U, and the entropy S 
of an ensemble of particles in volume V can be found from the 
thermodynamic relations P - -(dF/dV)T, S - -(dF/dT)v, 
and U - F + TS . The second-order thermodynamic functions 
are derived by differentiating these first-order ones. The decom- 
position (fT4l i induces analogous decompositions of P, U, S , the 
heat capacity Cy = {dS /dlnT)v, and the logarithmic derivatives 
XT = (<91nP/(91nr)y and Xp = -{dlnP/dlnVh- Other second- 
order function s can be expressed throug h these ones by Maxwell 
relations (e.g.. lLandau & Lifshitzlll980h . 



is the Fermi-Dirac integral. The internal energy is 



4 knTV 



[/3/2(^e,T) -HT/5/2(i'e,T)] . 



(21) 



Since we use V and T as independent variables, we need to 
find jUe(y, T). This can be done either by inverting Eq. (fT9] l 
numerically, or from the analytical approximation given in 
IChabrier & PotekhinI (Il998l) . Then the second-order thermody- 
namic functions are obtained using relations of the type 



df{x.,T) 
dT 

df(x.,T) 
dV 



df 




dT 


I: 




«e / 




V \ 





dXe 



(dnJdT\^ 
J (dns/dxe)T' 

V (dnjdx.h' 



(22) 
(23) 



We us e analytical app r oxima tions for Iy(xe,T), based on 
the fits of iBlinnikov et al.l d 19961) and accurate typically to a 
few parts in lO'*, with maximum error < 0.2% at t < 100 
(Chabrier & Potekhin 1998). These approximations are given by 
different expressions in three ranges of Xe- below, within, and 
above the interval 0.6 < Xe < 14. In particular, a t large Xe 
the Sommerfeld expansion (e . g. . IChandrasekhad 1 957t iGirif alcol 
1 19731) yield£] 

/vOfe, T) = —J— (jf (/i) + ^ t2 JP)(/}) + ...], (24) 



3. EOS of nonmagnetized Coulomb plasmas 

3. 1 . Ideal part of the free energy 

The free energy of a gas of A^i - n\Y nonrelativistic classical 
ions is 



F« =MA:Br[ln(M?/M,pi„)-l]. 



(15) 



where Mspi,, is the spin multiplicity. Accordingly, t/?^' 

\N,k^T, F® = n-,k^T, C®, = \N,k^, and;r?!,d = X^L = ^- 
In the OCP, Eq. ( fTSb can be written in terms of the dimension- 
less plasma parameters (Sect. |2} as 



^(i) o 16 

— — = 31n77- -InE- - In - -InM, 
Nik^T '2 In 



spin 



The free energy of the electron gas is given by 



(16) 



(17) 



where A^e = n^-V is the number of electrons and fi^ is the electron 
chemical potential without the rest energy nif^c^. The pressure 
and the number density are functions of //e and T: 



p(e) _ 



8 kuT 



3V^ Al 
4 



/3/2Cre,T) -I- -/5/2(ye,T) 



[h/2(Xe^T) + Tlij20Ce,T)] ., 



where ;ife = f^e/kBT, r = T IT,, and 
Iy(Xe,T)= — r—r^x 

Jo eXp{x-Xe)+l 



(18) 
(19) 

(20) 



where jl = XeJ - yUe/'WeC^ is the electron chemical potential 
(without the rest energy) in the relativistic units. 



[xf - ln(x + y)]/2. 



(0) , 



x^/3 

^3- 



5^2(//) = x'f/4 - 2x'/3 + 1.25 if l^ifi). 



(25) 
(26) 
(27) 
(28) 



Here, we have introduced notations x = yjil{2 + p) and y = 

Vl + x^ - 1 -H/i. At strong degeneracy, fi ^ ep/wieC^ - I, x ^ x,, 
and y y,.. In Paper II we also described an alternative expan- 
sion in powers of t, which allows one to avoid numerical can- 
cellations of close terms at small fi (we switch to it at /i < 0.1). 

The discontinuities of the IBlinnikov et aP (11996 ) approxi- 
mations for lyixe, t) at Xe - 0.6 and Xe = 14 are typically a 
few parts in 10"* at t < 10^, but they may reach ^1% for the 
second derivatives. This accuracy is sufficient for many appli- 
cations. Nevertheless, the jumps may produce problems, e.g., 
when higher derivatives are evaluated numerically in a stellar 
evolution code. In our calculations of white-dwarf evolution (to 
be published elsewhere), we smoothly interpolate between the 
two analytical approximations for the adjacent intervals near the 
boundary at cost of a slight violation of the thermodynamic con- 
sistency in the interpolation regions (this version of the EOS 
code is now also available at our web site). 

If a higher accuracy is needed, one can nu merically calcu- 
late tables of IviXe, t) (e.g., iTimmes & Amettl[l999 ) and inter- 
polate in the m with an a lgorithm that preserves thermodynamic 
consistency (Tim mes"&^Swestvll200Ql) and is available at MESA 
(Paxton etal. 2011). 



The multiplier 1 / V2 r"^' was accidentally omitted in Paper II. 
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3.2. Nonideal contributions 

3.2.1. Electron and ion liquids 

The contribution to the free energy due to the electron-electron 
interactions was studied by many authors. For the reasons 
explained in Paper II, we adopt the fit to Fee derived by 
llchimaru et all (Il987h (see Appendix|A]i. 

The ion-ion interactions are described on the base of the OCP 
model. In the liquid regime, the numerical results obtained for 
the OCP of nonrelativistic pointlike charged particles in differ- 
ent intervals of the Coulomb coupling parameter from T = 
to r ~ 200 by different numerical and analytical methods are 
reproduced by a simple expression given in Appendix IB. II The 
accurate fit for classical OCP is supplemented by the Wigner- 
Kirkwood correction, which extends the applicability range of 
our approximations to lower temperatures T ~ Tp.ln spite of the 
significant progress in numerical ab initio modeling of quantum 
ion liquids, available results do not currently allow us to estab- 
lish an analytical extensi on to still lower temperatures T <K Tp 
(see lChabrier et al.1l2002i for references and discussion). 

3.2.2. Coulomb crystal 

At r < Tn,, where is the melting temperature, the ions in ther- 
modynamic equilibrium are arranged in the body-centered cu- 
bic (bcc) Co ulomb lattice. In the harmonic approximation (e.g., 
lKittel[T963h . the free energy of the lattice is 

^'lat = f/o + f/q + f'th, (29) 

where Uq - NiCoiZef- ja^ is the classical static-lattice energy, 
Co ~ -0.9 is the Madelung constant, 

t/q = ^M^JWpMl (30) 

accounts for zero-point quantum vibrations, ui - {(jJka)^\\l<^p ~ 
0.5 is the reduced first moment of phonon frequencies, 

Fth = m.k^T <ln[l - exp(-Sa;fa,/^Br)]>ph (31) 

is the thermal contribution, are phonon frequencies, and 
{. . .)ph denotes the averaging over phonon polarizations a and 
wave vectors k in the first Brillouin zone. Here we do not sepa- 
rate the classical-gas free energy, therefore Fiat replaces F^^l+F.^:^ 
in Eq. 

Beyond the harmonic-lattice approximation, the total re- 
duced free energy /lat = Fi-^t/NiksT can be written as 

/lat = CoT +l.5ui}] + /th + /ah. (32) 

Here, the first three terms correspond to the three terms in 
Eq. (l29T l. and /ah is the anharmonic correction. The most accu- 
rate va lues of the constants Co and ui were calculated bv lBaikol 
(•2000) (see Appendix [B^ll. For f± = F±/NikBT, we use the 
highly precise fit of Baiko et al. (2001) (Appendix IB. 2b . In the 
classical limit 77 <s 1, it reduces to /h - 3 In 77 - 2.49389 - 
1.5mi77 -I- if I2A, where the term with u\ cancels that in Eq. (l32T i 
and the last term represents the Wigner-Kirkwood quantum cor- 
rection /p^ [Eq. (IB.3I 1I, which is the same in the liquid and solid 
phases tPoUock & Hansen. 1973.) . In the opposite limit T <s: Tp, 
we have /h - -209.3323 77"^ (here th e constant is given for the 
bcc crystal; for other lattice types, see Bai ko et al ] l200Th . 

Anharmonic corrections for Coulomb lattices were studied 
by many authors in the limits 77 — » and 77 — > 00, but only a 
few numerical results of low precision are available at finite 77 



values (see Paper II for references and discussion). In Paper II 
we have constructed an analytical interpolation between these 
limits, which is applicable at arbitrary 77 and is consistent with 
the available numerical results within accuracy of the latter ones 
(Appendix IB. 21 ). It should be replaced by a more accurate func- 
tion in the future when accurate finite-temperature anharmonic 
quantum corrections become available. 

3.2.3. Electron polarization 

Electron polarization in Coulomb liquids was studied by per- 
turbation dGalam & Hansenll976HYak ovlev & Shalvbkov 1989|) 
and hypernetted-chain (HNC) te chniques ( Chabrier & Ashcrofd 
ll990l:IChabrier& PotekhinI 19981: Paper I). The results have been 
reproduced by an analytical expression (Appendix IC.U . which 
exactly recovers the Debye-Hiickel limit for the weakly coupled 
(F <K 1) electron-ion plasmas and the Thomas-Fermi limit for 
the strongly coupled (F » 1) Coulomb liquids at Z » 1. 

For classical ions, the simplest screening model consists 
in replacing the Coulomb potential by the Yukawa potential. 
Molecular-dynamics and path-integral Monte Carlo simulations 
of classical liquid a rid solid Yukawa systerns were performed in 
sever al works (e.g., iHamaguchi et al 1 11 9971: iMihtzer & Grahariil 
l2006i) . However, the Yukawa interaction reflects only the 
small-wa venumber asymptote of the electron dielectric function 
(Ijancovici 1962; Galam & Hansen 1976). The first-order per- 
turbation approximation for the dynamical matrix of a classical 
Coulomb solid with the polarization corrections was developed 
by Pollock & Hansen (1973). The phonon spectrum in such a 
quantum crystal has been calculated only in the harmonic ap- 
proximation (Baiko, 2002,) , which has a restricted applicability to 
this problem (for example, it is obviously incapable to reproduce 
the polarization contribution to the heat capacity in the classical 
limit 77 ^ 0, where it gives Cy - 3M^b independent of the 
polarization). 

In Paper I we calculated Fj^ using a semiclassical perturba- 
tion theory of iGalam & HansenI d 19761) with a model structure 
factor, and fit the results by an analytical function of Xi- and 77. 
In Paper II we improved the 7/-dependence of this function so as 
to completely eliminate the screening contribution in the strong 
quantum limit 7/ «: 1 , because the employed model of the struc- 
ture factor failed at 77 < 1. The latter approximation is repro- 
duced in Appendix IC. 21 It can be improved in the future, when 
the polarization corrections for the quantum Coulomb crystal at 
77 < 1 will be accurately evaluated. 

3.2.4. Ion mixtures 

In Sects. I3.2.1I - [3".2.3I we have considered plasmas containing 
identical ions. In the case where several ion types are present in 
a strongly coupled Coulomb plasma, a common approximation 
is the "linear mixing rule" (LMR), 

/ex''*2-'^-^-(^>-'>"^^' ^^^^ 

where xj are the number fractions of ions with charge num- 
bers Zj, and Fj = ^eZf\ In Eq. (O, /ex = F^^/NiksT is 
the reduced nonideal part of the free energy, Fex - Fa for 
the "rigid" charge-neutralizing electron background and Fex = 
^^ii + ^ie + ^^ee for the polarizable background. The high ac- 
curacy of Eq. ( l33T l for binary ionic mixtures in the rigid back- 
ground was fi rst demonstrated by calculations in the HNC ap- 
proximation (iHansen et al.l 119771: iBrami et"aD 1 19791) and con- 
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firmed later by Monte Carlo simulations (iDeWitt et al.lll996t 
iRosenfeldlfToM iDeWitt & Slattervll2003h . The validity of the 
LMR in the case of an ionic mixture immersed in a polarizable 
finite-temperature e lectron background has been examined by 
[Hansen et al.l (Il977h in the first-order thermodynamic perturba- 
tion approximation and by lChabrier & Ashcroftl(ll990l) by solv- 
ing the HNC equations with effective screened potentials. These 
authors found that the LMR remains accurate when the electron 
response is taken into account in the inter-ionic potential, as long 
as the Coulomb coupling is strong (Fj > 1, Vj). 

However, the LMR is not exact, and Eq. ( l33T l should 
be replaced by the Debye-Hiickel formula in the limit of 
weak coupling (Tj 1, Vy). Even in the strong-coupling 
regime, the small deviations from the LMR are impor- 
tant f or establishing phase equilibria (see iMedin & Cummiiigl 
20101). The deviations fro m the LMR wer e studied by 



Brami et al. (1979); Chabrier & Ashcroft (1990); DeWitt et al.' 
(1996.) ; DeWitt & Slattery (2003); Potekhin et al.^ (2009a.b) for 
strongly coupled Coulorn b liquids and bv lQgata et al .1 (119931) and 
IDeWitt & Slattervl (12003') for Coulomb solids. 

The analytical expression that describes deviations from the 
LMR, A/ = / - f]M , in Coulomb liquids f or arbitrary coupling 
parameters Fj reads (iPotekhin et al.ll2009bh 



A/liq = 



V3 (i + fl<r)«)(i+/7<r>«) 



(34) 



where <Z*) = xjZ''., (T) = r^{Z^'^}, 6 is defined either as 

^'^ {zy/Hz^/2) 

for rigid electron background model, or as 

^ {Z(Z+l)^'^} «z2) + <Z))3/2 



(35) 



<ZV2) 



<Z)l/2 <Z5/2) 



(36) 



for polarizable background, and parameters a,b,a, and /3 depend 
on the plasma composition as follows: 



2.6 6+ 146^ 



l-a 



a = 



fe = 0.0117|^)fl, — -1. 



<Z2)l/5' 

3_ 

2^ 



(37) 
(38) 



For Coulomb solids, one should distinguish regular crys- 
tals containing different ion types and disordered solid mix- 
tures, where d ifferent ions are randomly distributed in regular 
lattice sites (Ogata et al. 1993). Each regular lattice type corre- 
sponds to a fixed composition, whereas random lattices allow 
variable fractions of different ion types. The free energy correc- 
tion A/ mainly arises from the difference in the Madelung en- 
ergies. It is generally larger for regular cryst als than fo r "ran- 
dom" crystals with the same composition. O gata et aP (Il993h 
performed Monte Carlo simulations of solid ionic mixtures and 
fitted the calculated deviation, A/soi, from linear- mixing predic- 
tion for the reduced free energy in a random binary ion crystal. 
[Medin & Cumming (2010) and Hughto et al. (2012) used this fit 
to study the phase separation and solidification of ion mixtures 
in the interiors of w hite dwarfs. We note, however, that the fit of 
iQgata et al.l d 1993b exhibits nonphysical features: for example, it 
is nonmonotonic as a function of Rz = Z2 /Zi at a fixed number 
fraction xt - \ - x\ for a binary ion mixture with Z2/Z1 > 2 and 
X2 < 0.5. A much simpler fit, which does not exhibit u nphysi- 
cal behavior, was suggested bv lDeWitt & Slattervl (l2003h . It can 



be written as A/soi - 0.00326 xiX2^2^^(r). However, the latter 
fit is valid only for relatively small charge ratios Rz < 3/2. We 
replace it by the expression 



A/sol ^ X1X2T1 Ag(x2,Rz), 
where 



(39) 



Ag(^2,^z) = 0.012 



-"^^ (1-R-^^)(1-X2+X2R'J') (40) 



-^2(1 - X2) 



and X — X2/Rz + (1 - Rz^) ■ Approximation (|40] | reason- 
ably ^wellre£roduces_both the results of Ogata et al. ( 1993 ) and 
De Witt & Slattervl (l2003l) for random two-c omponent ionic bcc 
lattices. For a multicomponent ion crystal, iMedin & Cummiiigl 
(I2OIOI) proposed the following extrapolation from the two- 
component plasma case: 



(41) 



A/,..^j;.,.x,r,A,(-^,|), 

where the indices are arranged so that Zj < Zj+i. 

4. EOS of a fully ionized magnetized gas 

4.7. Ions 



We consider only nondegenerate and nonrelativistic ions (for a 
discussion of the EOS of degenerate nuclear matter in stron 
magnetic fields see, e. g.. .Brodeiick et al. 20 00: Suh & Mathew 
120011) . In this case (cf. lPotekhin et alJI 19991) 



-(i) / -,9 



+ ln(l -e"^')- 



2 A^i^B^ 



(42) 



The last term arises from the energy of the magnetic moments of 
the ions in the magnetic field. 



f'spin = -A^i^B^ln 



sinh(gi^iMspin/4) 



sinh(gi^,/4) 



(43) 



where Mjpin is the ion spin multiplicity, and gi is the ^-factor 
(^i - 5.5857 for protons). For ions with spin one-half (Mspin = 
2), the expression in the square brackets in Eq. (l4Jt simplifies to 
[2 cosh(gi ^i/4)]. For zero-spin ions, such as "^He, ^^Fe, and other 
even-even nuclei in the ground state, Fspin = 0. 

The ion pressure obeys the nonmagnetic ideal-gas relation 
= n^k^T, but expressions for the internal energy and heat 
capacity are different: 



U 



(i) 



1 



A^i^B^ 2 ^ 



ef. - 1 2 



C 



(i) 
V,id 



A^i^B 



Here, the terms Mspin and Cspin arise from f spin. 



^.^i/4 



§i^iMspi„/4 



tanh(gi^i/4) tanh(gi^iMspi„/4)' 



- spin 



gi fiMspin/4 



, sinh(gi ^, /4) / \ sinh(gi ^iMsp,„ /4) 
They simplify at Mspin - 2: 



^spin 



tanh 



g.^./4 



cosh(g,^i/4) 



(44) 
(45) 

(46) 
(47) 

(48) 
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4.2. Electrons 
4.2.1. General case 

Thermodynamic functions of the electron gas in a magnetic 
field are easily derived from first principles ("Landau & Lifshita 
[l980). The number of quantum states per longitudinal momen- 
tum interval Ap, for an electron with given B-projections of the 
spin and the orbital moment and with a fixed Landau number » 
in a volume V equals VAp-/(47r'^alfi) jLandau & LifshitJ 19771) . 
Thus one can express the electron number density rig and the 
grand thermodynamic potential D.^f 



where pf^ is the zero-field Fermi momentum at the given den- 
sity. Equation i53[ can be written as pp - yii^cxb, where 



-P*'V as 

id 



An 



™ «=() o- 



dp. 



-oo exp[(e„(/7j) -jUe)/fcB7')] + 1 

yUe - eniPz) 



(49) 



In 1 -H exp 



dpz, (50) 



where Snip,) is given by Eq. (|7|, and 2^. denotes the sum over 
spin projections, which amounts to the factor 2 for n > 1 since 
we neglect the anomalous magnetic moment of electrons. This 
derivation equally holds in the relativistic and nonrelativistic the- 
ories. Equations ( |49] | and ( fSOl l can be rewritten, using integration 
by parts, as 

L^yy(i+2M'^^ '''^f"'^"\ (51) 



n=0 o- 



(52) 



where t„ = t/ V 1 + 2bn and Xn — Xe'^'^ Ths free energy 

pf^ is given by equations ( fTTI i. ( BTT i. and ( |52] ). 

The calculation of P\^, and their derivatives at given Xs 
and T can be performed using Eqs. (ISTt . (l52t and the same an- 
alytical approximations to the Fermi-Dirac integrals as for the 
nonmagnetized electron gas. The reduced electron chemical po- 
tential x<t at constant and T is found by numerical inversion 
of Eq. ( BTT i. Then the derivatives over T at constant V and over 
V at constant T are given by Eqs. (l22l i and (l23T l. We use this ap- 
proach in the current research, but we should note that for quan- 
tizing magnetic fields it is less precise than at Z? = 0. As men- 
tioned in Sect. 13.11 the inaccuracy of the employed approxima- 
tions for lyixe, t) is within a fraction of percent, but it grows for 
the derivatives. Since the first derivatives are already employed 
in Eq. ( BTT i. evaluation of the second-order thermodynamic func- 
tions such as xt or Cv involves third derivatives. Therefore, the 
error in the evaluation of these functions may rise to several per- 
cent. Such accuracy may be sufficient for astrophysical applica- 
tions, but otherwise one should resort to a thermodynamically 
consistent interpolation in numerical tables of the Fermi-Dirac 
integrals (Timmes & Arnett 1999; Timmes & Swestv 2000). 

Equations (ISTT l and ( |52] | can be simplified in several limiting 
cases considered below. 



4.2.2. Strongly quantizing and nonquantizing magnetic fields 

The field is called strongly quantizing, if most of the electrons 
reside on the ground Landau level. Then the electron Fermi mo- 
mentum equals 



p-e = iTi^alfin^ = (3/2)'''3 (flm/fle)^ /?f 



\2 



(53) 



xb 



2x] 
'3b 



30.2 



Z P6_ 



(54) 



is the relativity parameter modified by the field, and Xi- = 
/7pV'«eC (Sect. l2Tt . With increasing at constant B and zero 
temperature, the electrons start to populate the first excited 
Landau level when reaches hb = (7r^V2a^)"'. Therefore, 
the field is strongly quantizing at T « Tcyd and p < ps, where 
Tcyci = fiLoJkB ~ 1.3434 X IQ^Bn K and 



Pb = miriB/Z * 7045 (A/Z) B,'^ g cm" 



(55) 



The condition < ns can be written as a^ja^ < yfliZn) 
Then Eq. ( |53] | shows that in a strongly quantizing field /jp < 



/?p', except for densities close to the threshold n^. Thus 



is reduced, compared to its nonmagnetic value Tf\ by factor 



rp/r*"' = [(1 + xj^y^ - l]/iy, - l). in the nonrelativistic limit. 



,2^1/2 



Tp/Tf^ - ipp/Pp'Y, and the parameter - T/Tp turns into 



0B = 8yyj(9n^r,) ^ 0.166 Ooyirl 



(56) 



where 6*0 is the nonmagnetic value given by Eq. (|2]l. 

The opposite case of a nonquantizing magnetic field occurs 
at r » Tb, where Tb is the temperature at which the thermal ki- 
netic energy of the electrons becomes sufficient to smear their 
distribution over many Landau levels. It can be estimated as 
Tb = T'cyci, if P < Pb, and Tb = Tcyci/Tr, if p > Pb (a more 
sophisticated but qualitatively similar definition of Tb was intro- 
duced by Lai 2001). Then we can approximately replace the sum 
over Landau level numbers n by the integral over a continuous 
variable n. Integrating over n by parts, we can reduce Eq. (|52] | to 
Eq. ( fTSl ) and Eq. ( fSTT i to Eq. ( fT9b . i.e., to recover the zero-field 
thermodynamics. At p » p^, the electrons also fill many Landau 
levels, and the magnetic field can be approximately treated as 
nonquantizing. 

In the intermediate region, where the magnetic field is nei- 
ther strongly quantizing nor nonquantizing, the summation over 
n manifests itself in quantum oscillations of the thermody- 
namic functions with changing B and/or p, similar to the de 
Haas - van Alphen osci llations of magnetic susceptibility (e.g., 
iLandau & Lifshitjll980 ). The oscillations are smoothed by the 
thermal broadening of the Fermi distribution function and by the 
quantum broadening of the Landau levels (particularly, owing 
to electron collisions; see Yakovlev & Kaminker 1994., for ref- 
erences). Some examples of such oscillations will be given in 
Sect. El 

Figure[T]presents the p-T diagram of outer neutron-star en- 
velopes at two magnetic field strengths, B = 10'^ G and 10'^ G, 
assuming fully ionized iron (this assumption may be crude in 
the lower left part of the diagram). In the strongly-quantizing 
magnetic-field domain, bounded by pB and T^yd, the dependence 
7'p(p) is steeper than at B = 0, in agreement with Eq. (l56T l. 
The line Tj„(p) separates Coulomb liquid from Coulomb crys- 
tal. Near the lower right corner of the figure, where T <K Tp, 
the quantum effects on the ions become important (i.e., the ions 
cannot be treated as classical pointlike paiticles). In the lower 
left corner, at p < ps, the plasma is unstable to the phase separa- 
tion into the gaseous and condensed phases (this phase transition 
will be discussed in Sect. l7.3l l. 
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and the thermal correction to the free energy and internal energy 
AF = -AU = A^e Ae - VAP. (62) 

The leading contribution to the heat capacity is Cy^ - 2AU/T. 

As in the nonmagnetic case, Cy is propoitional to T at T — » 0, 
but with a different propoitionality coefficient. 

4.2.4. Strongly degenerate electrons in a strongly quantizing 
magnetic field 

If the magnetic field is strongly quantizing and the electrons 
are strongly degenerate (which coiTesponds to the triangular do- 
mains in Fig. [T] defined by conditions p < Pb and T < Tp), then 

fI^^ = [(1 + 4)'^' - l] N,m,c^ - (63) 
(2ny h(l+4)'^'-ln(-^B + (l+4)'^')l- (64) 



' 1 ' 

- 






/ I 


1 

- 
- 










- 






" / 










^ / 






^ cycl \^ 12 ^ / 




<s/ 




- 




, 


" r 


\ 








= 1000)' 


O 1 

o 

II ' 




.J 


! \ 
' I a" 

1 1 , 1 < 1 , 


CM ' 

s 1 

a 1 
1 , 


<?! 
, 1 1 


. 1 



5 6 7 8 
log p (g cm-3) 



10 



p(e) _ 



Fig. 1. Characteristic density-temperature domains at B = 
10^^ G (blue online) and 10'^ G (red online) for fully ionized 
iron. Solid lines indicate the Fermi temperature as function of 
density, the dotted line shows the plasma temperature, the dot- 
dashed line shows the melting temperature as function of den- 
sity, short and long dashes delimit the domains of strongly quan- 
tizing magnetic field and of magnetic condensation, respectively, 
and the heavy dots mark the critical point for the condensation 
(Sect.ESll. 

4.2.3. Strongly degenerate electrons 

If the electrons are strongly degenerate, then one can apply the 
Sommerfeld expansion 

(Sect.Ol and obtain Ff^ k. Ff + AF, 

(q) (q) 

where F^ — epA^e - ^ zero-temperature value and 

AF is a thermal correction. According to Eqs. (l24t and (l52t . the 
zero-temperature pressure Pq is 



^o'=^ZZ(i+2MJr;#„), (57) 

«=() IT 

where P^ = m\c^ IfP' - 1.4218 x 10^^ dyn cm"^ is the relativis- 
tic unit of pressure, Wniax is the maximum integer n for which 
p\{e) > 0, and e„ = ep/mgC^ -H 1 - Vl + 2bn. According to 
Eqs. (l24l i and ( fSTT l. the Fermi energy ep is determined by the 
condition 

= (^F 2^ S 2M'%i(f„). (58) 

n=() iT 

In order to obtain the chemical potential yUg = ep + Ae with frac- 
tional accuracy ~ Xe^^ we retain two terms in Eq. (l24l i. substitute 
it into Eq. ( BTI ). approximate I^"\ji„) in the vicinity of //„ = e„ 
by 

I["\fi„) ^ P^\~e„) + Pr'\e„) Al, (59) 

where /}„ = XnT and Ae = Aejm^c^, and drop the higher-order 
terms containing (r^Ae). Then 



Ae ~ - 



k't' 2^(1 + 2M"'^'^f/2(^") 



6 2:;;:o2-(i + 2M'^'0^«) ' 

The thermal correction to the pressure equals 



(60) 



AP = J ^(1 + Ihn) I^^lien) + ^ l[%ien)] , (61) 



27r2 



In the nonrelativistic (xfi <k 1) and ultrarelativistic (x^ » 1) lim- 
its, we have Pf ^ P, b x\l6Tr- oc nl and P'^^^ ^ P, b x\IA7t^ cc nl, 
respectively. Compared with the nonmagnetic case (Papers I and 
II), the dependence *(«e) is steeper, but P is lower everywhere 
except for «e ~ ng- Thus, a strongly quantizing magnetic field 
softens the EOS of degenerate electrons. 

The thermal coiTections (l60t - (l62t simplify to 



Ae^ - 



6(1+4)1/2x2' 



AP = P, 



bT^ 



2 + xi 



12 (l+x2)i/2x/ 



AU AF „ br^ (1 +x|)i/2 



V V 12 ' 



Xb 



C 



(e) 



;r2T(l+x2)i/2 



(65) 



(66) 



(67) 



(68) 



The last equation differs from to the nonmagnetic case (Paper II) 
by the replacements of by xb and by 7r^/3. 

4.2.5. Nonrelativistic limit 

In the nonrelativistic limit (pp <K mgC and T <K T^), Eqs. ( BTT i 
and ( |52] | simplify to 

'Ti ^ n.cr 



27r3/2fl^i, 



In the nondegenerate regime (T » Jp), one has /vOf) ~ S''^ ^(v + 
1), were r(v -H 1) is the gamma-function. Then Eq. ( |69l ) yields 
P^a = ne^B^ and 

;re = ln(«e^e/2) " ln(f /4) + ln[tanh(f/4)] (70) 

which provides the free energy P^* - (Xe ~ 1) A^e^B^. 

As follows from Eq. (iTOl i. the reduced internal energy 
U^JNiksT and heat capacity Cy-^JNik^ of the Boltzmann gas 
decrease with increasing ^. In a strongly quantizing magnetic 
field (^ » 1), they tend to 1/2 instead of 3/2, because the gas 
becomes effectively one-dimensional, with the only kinetic de- 
gree of freedom along the magnetic field. 

In the nonquantizing field (^ <s; 1), the two last terms in 
Eq. ( iTOb cancel out, so that the standard expression F^^ - 
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NeksT [ln{nsAl/2) - 1] is recovered. In the strongly quantizing, 
nondegenerate regime (p < ps and Tp <K T <K Tcyd), the last 
term of Eq. (iTOl i vanishes, which yields 



Flf ^ N,kjiT[ln(27TaiA,n,) - l]. 



4.3. Thermodynamic and kinetic pressures 



(71) 



The above expressions for pressure of a magnetized gas of 
charged particles are based on the principles of thermody- 
namics (iLandau & Lifshitd Il980h . according to which P - 
-{dF / dV)T,B- Alternatively, the pressure can be calculated from 
the microscopic dynamics as the sum of the changes of ki- 
netic momenta of all particles reflected off a unit surface per 
unit time. The result of the latter calculation (the "kinetic pres- 
sur e" p^ ") depends on the orientation of the surface relative to 
B (ICanu to & Chiu 1968). If the surface is perpendicular to B, 
then one gets the kinetic pressure P^" = P, which acts along the 
field lines (the "longitudinal pressure"). If the surface is parallel 
to the field, one gets a d ififerent ("transverse") kinetic pressure, 
which can be expressed dBlandford & Hernquist|[T982b as 



P^i" = -Q/y + B (dn/dB)v,T ^P- MB, 



(72) 



where M — -dQ/dB is the magnetization. 

In order to resolve the apparent paradox, one should take the 
magnetization current density jm = c V x M into account; in 
case of boundaries, this volume current should be supplemented 
by the s urface current cM x B/B (see, e .g., iGriflithlll999l) . As ar- 
gued bv iBlandford & Hernquisl (Il982h . if we compress the elec- 
tron gas perpendicular to B then we must do work against the 
Lorentz force density j,n x B/c, which gives an additional con- 
tribution to the total transverse pressure and makes it equal to P. 
Because this point still causes confusion in some publications, 
let us illustrate it with a simple example. 

If the pressure were anisotropic, then one might expect an 
anisotropic density gradient in a strongly magnetized star. Let 
us consider a small volume element in the star, assuming that 
we can treat B, T, and gravitational acceleration g as constants 
within this volume, and z axis is directed along g. Hydrostatic 
equilibrium implies that the density of gravitational force, pg, 
is balanced by the density of forces created by plasma particles. 
The crucial point is that the magnetization contributes to this 
balance. 

Let us compare the cases where B is parallel and perpendic- 
ular to g. In the first case, the z-component of Lorentz force is 
absent, and we get the standard equation of hydrostatic equilib- 
rium: pg = dPjj^'"/dz = dP/dz. In the second case, the gradients 

of the kinetic pressure dP^"/dz and of the Lorentz force density 
BdM/dz act in parallel. In the constant and uniform magnetic 
field, dM/dz is not zero, but is related to the density gradient: 



dM 
dF 



dMjp, T, B) dp 
dp dz 



d^Qip, T, B) dp 
VdpdB dF' 



(73) 



Then the equilibrium condition takes the form pg - dP^'^/dz + 
BdM/dz = dP/dz, which is the same as in the first case. 
Furthermore, one can express P through p using some EOS. In 
the considered example, dP/dz - [dP(p, T, B)/dp] dp/dz, so that 
dp/dz is also the same in both cases. Thus, the stellar hydrostatic 
profile is determined by the isotropic thermodynamic pressure P, 
which automatically includes magnetization. 



5. Magnetic effects on the EOS of a Coulomb liquid 

5.1. Electron excliange and correlation 

The eff'ects of a magnetic field on the contribution to the free 
energy due to electron exchange and correlation were studied 
either in the regime of strong degeneracy and strongly quantiz- 
ing magnetic fields dPanz & Glas ser"1971' 'Baneriee etal][T974l: 
iFushiki et al.]|1989l: see also lMorbec & Capelle 20081 for an in- 
structive discussion of the previous results and the inclusion 
of the second Landau level contribution), or at low densities 
( Alastuey & Jancovici 1980; Cornu 1998; Steinberg et al. 199^, 
2000). In a previous work (iPotekhin et al.lll999i) we suggested a 
modification of the field-free expression for Fee, which matches 
available exact limiting expressions, including the cases of non- 
quantizing, strongly quantizing degenerate, and strongly quan- 
tizing nondegenerate plasmas. The modification consists in re- 
placing Ffie{0, Fe) by Fee(0*, Fg), whcre 



+ 0b 



Ob , cosh(^/2) tanh(^/4) arctanh^ 

i -H — exp(-ty£ ) 



, (74) 



[cosh(^/4)]2 ^/4 



^ = [1 - (4/^) tanh(^/4)] ^1^, and 6*0 and 9b are given by Eqs. ^ 
and ( l56] i. respectively, at fixed iie and T. 

5.2. Wigner-Kirkwood term 

For the same reason as in Sect. 13.2.11 the treatment of the 
quantum effects in the ion liquid is restricted by the Wigner- 
Kirkwood term. Its expression in an arbitrary magnetic field was 
obtained bv Alastuev & Jancovici (1980,) : 



1 



<2) 



24 



4 8 1 

tanh(^i/2) " ^ ^ 3 



(75) 



The function in the square brackets monotonously varies from 1 
at — » to 1/3 at oo, reflecting the effective reduction of 
the degrees of freedom of the classical ion motion from d - 3 
at B = to d = 1 for a strongly quantizing field. At small ^i. 



/f * (7,^/24) (l-ff/90). 



5.3. Electron-ion correlations 

Using the linear re sponse theory in the Thomas-Fermi limit, 
iFushiki et al.l (Il989l) evaluated the electron polarization energy 
for a dense plasma in a strongly quantizing magnetic field at zero 
temperature, assuming that the ions remain classical (unaffected 
by the field). A comparison with the analogous zero-field result 
shows that the strongly quantizing magnetic field (ymr? > 2.23) 
increases the polarizati on energy at high de nsities (r^ <K 1) by a 
factor of 0.8846 ylrf (Potekhin et al." 1999'). 

Recently, Sharma & Reddv (2011) calculated the screening 
of the ion-ion potential due to electrons in a large magnetic field 
B at r = 0, using the one-loop representation of the polarization 
function. Their results for the strongly quantizing magnetic field 
show that the screening is anisotropic, and the screened ion po- 
tential exhibits Friedel oscillations with period nh/ pp in a cylin- 
der of a radius ~ nh/pp alon g the magnetic field line that passes 
through the Coulomb center. ISharma & Reddvl suggest that this 
long-range oscillatory behavior can affect the ion lattice struc- 
ture. However, finite temperature should damp these oscillations, 
so that they are pronounced only at T <sc Tp, i.e., deep within the 
triangular domains formed by the lines Tp and p^ in Fig. [1] At 
the typical pulsar magnetic fields B ~ 10'^ G, this requires an 
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unusually low temperature of the neutron-star crust. On the other 
hand, the conditions T <K Tp and p < ps can be easily fulfilled 
in the outer crust of magnetars at B ~ 10'^ G (cf. Fig. [T] and the 
top panel of Fig. [8]l, but in this case the Friedel oscillations are 
strongly suppressed because the electrons are ultrarelativistic. 

To the best of our knowledge, the magnetic eff'ects on the 
electron polarization energy have not been calculated at finite 
temperatures or in the case where the field is not strongly quan- 
tizing. In view of the limited scope and limited applicability of 
the available results on the magnetic effects, we use the nonmag- 
netic expression for Fi^ in our code (AppendixICli. 



6. Harmonic Coulomb crystals in the magnetic field 

The magnetic effects on Coulomb crystals have been studied 
only i n the harmonic approximation. iNa^ai & Fu kuvamal( ll982i 
Il983h calculated phonon spectra of body-centered cubic (bcc), 
face-centered cubic (fee), and hexagonal closely-packed (hep) 
OCP lattices. They compared the energies of zero-point vibra- 
tions at different values of parameters /3 and a nd found con- 
ditions of stability of every lattice type. However. Baikd (l2000i 
[2009) noticed that their choice of the magnetic-field direction 
did not provide the minimum of the total energy. 

lUsov et aP (|T980) obtained the equations for oscillation 
modes of a harmonic OCP crystal and studied its phonon spec- 
trum in a quantizing magnetic field in several limiting cases. 
These authors discovered a "soft" phonon mode with dispersion 
relation ojka ^ near the center of the Brillouin zone, which 
leads to the unusual dependence of the heat capacity of the lat- 
tice Cv_\M °^ T^^^ r — » instead of the Debye law Cvm T^. 
IUsovetal.l(IT980l) argued that a strong magnetic field should in- 
cre ase stability of the ciy stal. 

iBaikol (120001 |2009|) studied the magnetic effects on the 
phonon spectrum of the harmonic Coulomb crystals and calcu- 
lated its energy, entropy, and heat capacity. We have found that 
his results can be approximately reproduced by analytical ex- 
pressions presented below. 



6. 1 . Thermal phonon contributions 

Without magnetic field, the thermal phonon contribution /th to 
the reduced free energy of a Coulomb crystal F/NiksT is a func- 
tion of a single argument described by a simple analytical ex- 
pression (Baiko et al. 2001). The magnetic field introduces the 
second independent dimensionless argument The functional 
dep endence of th ermodynamic functions on i] and /3 is not sim- 
ple. [Ba&^'llQQ^ identified five characteristic sectors of the r]-/3 
plane: 

1. ?7 < 1 and /3 < rj^^ - weakly magnetized classical crystal, 
T] > I and p <ri^^ - weakly magnetized quantum crystal, 
77 < 1 and y6 > 77"' - strongly magnetized classical crystal, 
T] > /3 > Tf^ - strongly magnetized quantum crystal, 
/? > 77 > 1 - very strongly magnetized quantum crystal. 
Note that the condition p > rf^ \s equivalent to > 1. Thus, the 
magnetic field strongly affects the thermodynamic functions of 
a Coulomb crystal when ticjci > k^T, that is the same condition 
as for the gas and liquid phases. 

For astrophysical applications, we have constructed an an- 
alytical representation of the EOS of the magnetized Coulomb 
crystal, which is asymptotically exact in every of the five sectors 
far from the ir boundaries, exactly recovers the nonmagnetic fit of 
iBaiko et al] ([2001) in the limit y6 — > 0, and reaches a reasonable 
compromise between simplicity and accuracy. 



2. 
3. 
4. 
5. 




-2 -1 
log T-Ap 

Fig. 2. Thermal phonon contribution to the reduced internal en- 
ergy Mth = U^h/Nik^T as a function of log(r/7'p) = -log 77 at 
yS = noj.i/kBTp = 0, 1, 10, 100, and 10^ (numbers near the lines). 
Analytical approximations (l76T l (dotted lines) and (fTST l (short- 
dashed lines) are compared with the numerical results of iBaikol 
(I2009h (sohd lines for/? = 1, 10, and 100). 
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Fig. 3. Thermal phonon contribution to the reduced entropy 
■Sth = Sth/NiksT as a function of log(r/rp) at yS = hLOd/ksTp = 
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Fig. 4. Thermal phonon contribution to the reduced heat capac- 
ity Ci/iat/M^B^ as a function of log(r/rp) at yS = HdJajk^T^ - 0, 
1, 10, 100, and 10^ (numbers near the lines). Analytical approx- 
imation ( l80l l ( short-dashed lines) is compared with the numer- 
ical results of lBaikd (l2009l) (solid lines for yS = 0, 1, 10, and 
100). The dotted lines correspond to the first term on the rh.s. of 
Eq. dlOll. 



The term f± in Eq. ( |32] | can be rewritten as /th - u± - .Sth, 
where u± - Uth/NiksT and ith = Sth/Niks are the thermal con- 
tributions to the reduced internal energy and entropy. We approx- 
imately represent Mth by the function 



.,(0) 



M - 



-2/5 



lA / yi + 24/1]^ 

I + 15/3/1] + 



1+0.5 [1 +(3/^i)5/2] 
where i// = 12.5 (/3/t])^'^ + 1 19 (J3/T]f, = /3t], 

and ith by the function 



(76) 



4°)+ in 



1 + 



20.8 i/S/T])^'^ + 122.5 (/3/T]f 
(l+7.9/?/77)(l+ 42/772) (l+A-i)4j 



(77) 



In these equations, i/?' and ijf'' are the values of Mth and .Sth at 
j8 = 0. Equations ( f76] l and dTTl ) exactly reproduce the known 
asymptotic limits: Mth = 3 in the classical nonmagnetic limit 
(77 <K y6 <K 1), Mth = 2 in the classical magnetic limit {rj 

' <K 1), Mth = 1 in the case where /? » 77 » 1, and 
Mth = 0.6sth ^ 1]^^^^, if 77 —> 00 at y6 = constant. 

The functions u{ri,/i) and s{ri,/i) are displayed in Figs.|2]and 
[3] Their accuracy is seen from a comparison with the numeri- 
cal results ([Baiko 2009), also shown in the figures. However, if 
the complete consistency of different thermodynamic functions 
is required, Eqs. (l76t and (l77t should not be used directly, but 
should be first combined into /th - u - s. Then one can cal- 
culate thermodynamic functions by differentiating the function 
/th(77,y6). In this way we obtain, for example. 



Mth - 



^th 



— u + Am, 



(78) 



Cy,th / ds 



s + Am, 



NikB 

where Am = 



5 In 77^^ 
ds 



dAu 
5 In 77 

du 



d\nri/„ \5in77 



(79) 
(80) 

(81) 



Note that the relation between the phonon contributions to pres- 
sure and internal energy, 2f th^ = t^th, which is standard for a 
nonmagnetized harmonic crystal, is invalid in the strongly mag- 
netized crystal, because both dimensionless arguments 77 and /3 
depend on density. 

Approximations (iTSTi - (I8OI1 are shown in Figs. |2]-|4] Their 
reasonable behavior beyond the range of available numerical 
data is demonstrated by plotting them also at larger ^6 = lO-' 
and 10^. 



6.2. Zero-point vibrations 

Because motions of the ions are confined by the magnetic field 
in the transverse direction, they exhibit quantum oscillatio ns in 
the ground state (lLandaulll93a iLandau & Lifshitdll977h . The 
energy of these oscillations is TiUa/l for every ion, which gives 
the term (\/2 in Eq. (l42l i. In a Coulomb crystal, the motion of an 
ion is confined in an effective potential well, centered at its equi- 
librium lattice site. The total energy of the zero-point quantum 
lattice oscillations is given by Eq. ( l30l l. 

In the case where the crystal is placed in a magnetic field, 
Uq includes contributions due to both, magnetic and lattice con- 
finements of the ion motion. However, since the magnetic con- 
tribution Nih(jL>a/2 is common in all phase states, we take it as 
the zero energy point in our code and separate it from the lat- 
tice contribution that is specific to the solid phase. Then Eq. ( l30t 
becomes 



3 , 1 



(82) 



and in Eq. (l32b we have 1.5mi77 = 1.5mj77 + ^i/2, where we have 
defined 



Mj - Ml -yS/3. 



(83) 



The reduced frequency moment u\ still depends on y6, because 
the character of ion vibrations is affected by the magnetic field 
(they become essentially one-dimensional if B is extremely 
large), but the latter dependence is relatively weak. Havi ng ex- 
tracted u \ from the available numerical results for ui (iBaikd 
2009; Ba iko & Yakovlevll2012h . we can represent it by the sim- 
ple interpolation 



«'i(y6) = 



1.27/3' 



.9/8,, 



1 + 1.27/3'^/8 



(84) 



where m" is the zero-field value and m™ is the limit of u'^ifi) at 
/3 — > 00. Only one of the three phonon branches contributes to u[ 
in the latter limit, therefore m^ ~ m*jV3. For the bcc crystal, m^ = 
0.51 13875, whereas u'^ varies between 0.18 and 0.19 depending 
on the orientation of the lattice in the magnetic field. In Fig. |5] 
we show M| an d the lo garithm (base 10) of mi versus /?. 

lUsov et al ] (11980') noticed that the energy of a crystal de- 
pends on its orientation in a strong magnetic field. However, nu- 
merical calculations (Baiko 2000, 2009) show that this depen- 
dence is very weak. For example, the difference Ami between 
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0.01 0.1 1 10 

Fig. 5. Reduced first moment of phonon frequencies m'j (solid 
line) and logwi (dashed line) as functions of fi - fiajci/Tj,. The 
dotted line shows the asymptote log(J3/3). 

the values of ui for two orientations, where the field lines con- 
nect an ion with its nearest neighbor in the first case and with a 
next-order nearest neighbor in the second case, is approximately 

«15/4 

" (1 +^3/2)5/2 (^"l)"-" (85) 

with saturation level (AwOmax = 0.0064 for the bcc lattice. 

6.3. Comparison with the Baiko-Yakovlev fit 

After the present work had been completed, we beca me aware of 
an independent study bv Baiko & YakovlevI (l2012h . who devel- 
oped another set of approximations for the free energy of the har- 
monic Coulomb crystal in a magnetic field. They presented the 
free energy as a sum of three terms corresponding to the contri- 
butions from each of the three phonon modes in the bcc crystal. 
Thus every of these terms has a clear physics meaning, while our 
fitting expressions give only the total contribution, which cannot 
be easily decomposed to three parts corresponding to the sepa- 
rate phonon modes. 

Unlike our fit, the fit of iBaiko & Yakov lev (20 121) does not 
exactly reproduce the very accurate results of B aiko et al.l(l2001h 
in the limit y6 — > 0. At finite p, both sets of fitting expressions ac- 
curately reproduce the asymptotes at T <K Tp and T » Tp and 
have similar accuracies within several percent in the intermedi- 
ate range O.lTp/yS < T < IQTp. Meanwhile, our approximation 
is simpler: the Baiko-Yakovlev approximation contains 27 inde- 
pendent numerical fitting parameters, whereas our fits (fTFi and 
(TTTl i contain together only 9 such parameters. 

7. Examples and discussion 

7.1. Thermodynamic functions 

Characteristic features of the EOS can be seen in Fig. |6] Here, 
we have chosen the plasma parameters that are typical for outer 
envelopes of isolated neutron stars: we consider fully ionized 



iron (Z = 26, A = 56) at T = 10^ K and B = 10'^ G (for il- 
lustration, the density range is extended to p < 10^ neglecting 
the bound states that can be important in this p-T domain). We 
plot the normalized pressure p - P/riiksT, entropy S /NikB, heat 
capacity cy = Cv/Niks, and logarithmic derivatives of pressure, 
Xp andxT as functions of density. Dashed lines show these func- 
tions in the absence of quantizing magnetic field. The vertical 
dotted lines marked by numbers separate different characteristic 
domains, consecutively entered with increasing density: onset 
of electron degeneracy at B = (Tp"* - T) and ?A B - lO'^ G 
(Tp = T), population of excited Landau levels (p - pg), melting 
point with formation of a classical Coulomb crystal {T^ = T), 
quantum effects in the crystal (Tp = T). 

At low densities, the ideal-gas values are approached: p = 
1 +Z,Xp =Xt = l,ci/ = (3 H- 3Z)/2 at B = 0, andc^ = (3-i-Z)/2 
at B = 10'^ G. The latter difference is caused by the fact that at 
p < Pb the electron gas is effectively one-dimensional due to the 
strong magnetic quantization. 

With increasing density, the reduced pressure p first de- 
creases below its ideal-gas value due to the Coulomb nonide- 
ality and then increases due to the electron degeneracy. The in- 
crease occurs earlier at zero field than in the strong magnetic 
field, because of the delayed onset of the degeneracy (Sec. |4]i. 
When p > Pb ~ 1.5 x 10^ g cm""*, the thermodynamic functions 
approach their zero-field values. The gradually decreasing oscil- 
lations correspond to consecutive filling of the electron Landau 
levels. The magnetic field B = 10'^ G does not affect the ion 
contributions in Fig. |6] because it is nonquantizing for the iron 
nuclei at r = 10^ K (^i = 0.00342). 

The liquid-solid phase transition occurs in Fig. |6] at p » 
8.25 X 10"^ g cm"-', where we adopt the classical OCP melting 
condition F - 175.2 (Paper I). With further increase of density 
(p > 10''), the degeneracy becomes so strong that the energy and 
pressure are nearly independent of T, wAxt strongly decreases. 
The normalized heat capacity gradually tends to its value cy -1> 
characteristic of the classical simple crystal. At still higher den- 
sity, the ion motions become quantized (Tp » T), which leads 
to the further decrease of the heat capacity and the entropy. 

7.2. Melting 

The electron polarization, ion quantum effects, and quantizing 
magnetic field can shift the melting temperature. The lower 
panel of Fig. Q shows the Coulomb coupling parameter F at the 
melting (that is, the value F„i at which the free energies of the 
two phases are equal to each other); the upper panel displays the 
difference between the internal energies in the liquid and solid 
phases at the melting point (the latent heat 2m - Ua^ - Usoi), di- 
vided by NiksT. We plot the data for fully ionized '^C and ^^Fe 
at B = and 10'^ G for carbon, B = 0, 10"* G, and lO''* G 
for iron. The density range shown in the figure is typical for the 
outer envelopes of a neutron star and is also relevant for white 
dwarfs. 

The position of the melting point is very sensitive to the ac- 
curacy of the free energies of the Coulomb liquid and crystal 
(see, e.g.. Paper 1). The polarization and quantum corrections to 
the classical OCP free energy are not known sufficiently well for 
finding the position of the melting point in the whole interval 
of densities shown in Fig. Q The dot-dashed curves in this fig- 
ure correspond to the domain where T < 0.5 Tp. Here, the per- 
turbation theory for the quantum effects in the liquid phase be- 
comes progressively inaccurate. The dashed curves correspond 
to the domain where the binding energy of the hydrogenlike 
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Fig. 7. Characteristics of the melting transition of nonideal car- 
bon and iron plasmas at different field strengths B (marked near 
the curves). Lower panel: the value of the Coulomb coupling 
parameter F at the melting point as function of mass density 
p. Upper panel: normalized latent heat per ion at the melting 
transition. The dot-dashed and dashed segments of the curves 
correspond to the domains of nonperturbative quantum effects 
(T < 0.5 T-p) and electron response (Z^ Ry > 0.1 ep), respec- 
tively. The dotted horizontal lines mark the OCP values. The 
filled and open circles mark the positions of the real and virtual 
condensed surfaces (see text in Sect. 17.3b . 



ion exceeds 10% of the Fermi energy, Z^Ry > 0.1 ep, which 
corresponds to > O.baojZ. Here, the perturbation treatment 
of the electron polarization in the crystal starts to be inaccu- 
rate. In addition, the position of the melting point cannot be 
traced to F < 100-120, because the available results for the 
anharmonic corrections to the free energy of a Coulomb crystal 
(Appendix IB. 2b are accurate only at larger F. Nevertheless, we 
can evaluate F^ in a certain interval of densities for each type of 
ions (the solid segments of the curves in the figure). 

The values of gm in Fig. |2]roughly (within a factor of two) 
agree with the OCP value Q^^ = Q.lAeNik^T at F,,, = 175.2 
and with the values used in theoretical models of white dwarf 
cooling (e.g., Hansen 2004 and references therein). Most of the 
neutron star cooling models currently ignore the release of the 
laten t heat at the c rystallization of the neutron star envelopes 
(e.g.. lYakovlev & Pethick 2004, and references therein). 

Figure Q shows that the strong magnetic fields tend to de- 
crease Fm and thus stabilize the Coulomb crystal, in quali- 
tative agreement with p revious co njectur es ( Rudermaiilll97ll : 
iKaplan & Gla"sse3IT972L lUsov et alJll980l lLaill2001h . At densi- 
ties p ~ 10^ - 10^ g cm~^ correspon ding to the "sen 
in the neutron-star cooling theory (lYakovlev & Pethici 
the stabilization proves to be significant at the magnetar field 
strengths B = lO'"*- 10'^ G. The results for ^'^Fe in this B inter- 
val, shown in the lower panel of Fig. [T] can be roughly (within 
10%) describedby the formula Fn,(B) « Fm(0)/(1 -1-0.2/3). At the 
typical pulsar field strengths, B - 10'^- 10'^ G, the effect is no- 
ticeable at lower densities. However, these conclusions remain 
preliminary in the absence of an evaluation of the magnetic-field 
effects on the anharmonic corrections. 

In view of the limited applicability and incompleteness of 
the evaluation of F^ with account of the quantum, polarization, 
and magnetic effects, in applications we use the classical OCP 
value Fm = 175.2 as the fiducial melting criterion. 

7.3. Magnetic condensation 

iRudermanI (Il97ll) suggested that the strong magnetic field may 
stabilize molecular chains (polymers) aligned with the magnetic 
field and eventually turn the surface of a neutron star into the 
metallic solid state. Later studies have provided support for this 
conjecture, although the critical temperature T^nt, below which 
such condensation occurs, remains very uncertain. Condensed 
surface density ps is usually estimated as 



it y strip " 
M l2004h . 



■ 561 ^AZ-^^-^B\f 



gem 



(86) 



where ^ ~ 1 is an unkno wn numerical factor, which ab sorbs 
the theoretical uncertainty (lLaill2001l: iMedin & Laill2006 ). The 
value ^ = 1 correspo nds to the EOS provided by the ion- 
sphere model (Salpete nll961l). which is close to the "uniform 
model" of Fushiki et al.1 (Il989l) . For comparison, the results of 
the zero-temperat ure Thomas-Ferm i model for ""^Fe at 10'° G < 
B < 10'^ G ( IRognvaldsson et al.lll99 3) can be approximated 
(within 4%) by ps,^ with ^ ^ 0.2 + 0.01 /fi' ,^'', whereas the 
finite- temperature Thomas-Fermi model of IThorolfsson et al.l 
(Il998h does not predict magnetic condensation at all. Our 
EOS f or partially ionized hydrogen plasmas in st rong magnetic 
fields (Potekhin et al. 1 999l:IPotekhin & Chabrier 2004) exhibits 
a phase transition with Tcm ~ 3x10^ ^12^'' ^ ^'^'^ critical density 
Peril ~ 143 B\-^ ^ g cm~^ ~ ps.o 25 at 1 < B p < 1 0^ According to 
another study (iLai & Salpetedll997l: lLaill200ll) . T^^it for hydro- 
gen is_severaltimessmal^ 

iMedin & Lail (l2006l) performed density-functional calcula- 
tions of the cohesive energy Qs of the condensed phases of H, 
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He, C, and Fe in strong magnetic fields. A comparison with pre- 
vious density-functional calculations of other authors prompts 
that Qs may vary within a factor of two at B12 > 1, depending on 
the approximations (see iMedin & Lai 2006 for r eferenc es and 
discussion). In a subsequent study, Medin & Lail (l2007l) calcu- 
lated the equilibrium densities of saturated vapors of He, C, and 
Fe atoms and polymers above the condensed surfaces, and ob- 
tained Tcrii at several values of B by equating the vapor density 
to ps. Unlike the previous authors, iMedin & Lail (l2006t 120071) 
have taken the electronic band structure of the condensed mat- 
ter self-consistently into account, but they did not allow for the 
atomic motion across the magnetic field and mostly neglected 
the contributions o f the excited a tomic and molecular states in 
the gaseous phase. iMedin & Lail obtained ps assuming that the 
linear molecular chains form a rectangular array with sides 2R 
in the plane perpendicular to B, and that the distance a be- 
tween the nuclei along B in the condensed matter remains the 
same as in the gaseous phas e, so that Ps = mJA aR^ (Medh] 
I2OI2I) . Using Tables III- V of iMedin & Lail ( l2006l) for a and R, 
we can describe their results for the surface density of '^C at 
1 < B12 < 1000 and ^^''Fe at 5 < Bn < 1000 by Eq. (|86]l with 
f = 0.517 -H 0.24/fiJ^^ + 0.0 11 and ^ = 0.55 + 0.1 1, respectively. 

iMedin & Lail (12007 ^) found that the ciitical temperature is 
2"crit ~ 0.082s/^B- Their numerical results for He, C, and Fe 
can be roughly (within a factor of 1.5) described as r^rit ~ 
5 X io4z i/4g3/4 ^ ^ g ^ ^QQQ P^j. comparison, the 

results of iLai & SalpeteJ (1 19971) for H at 10 < Bn < 500 sug- 
gest Tcrit ~ 1 .6 X 10^ Bjj' K. The discrepancies between different 
estimates of ps and Tait reflect the current theoretical unceitain- 
ties. 

The filled dots on the curves for magnetized carbon in both 
panels of Fig. |7] coiTespond to the condensed surface posi- 
tion in the fully ionized plasma model. In this model, T^rit ~ 
2.5 X 10^ Z"-' 5**2* K and pcrit ~ ps.().47. With decreasing temper- 
ature T below Tcrit, the surface density increases and tends to 
the hmit psj given by Eq. (l86l l as ps ~ Ps,i/[1 + 1-1 {T/Tcm)^] 
(cf. Fig. [TJ. At smaller densities, there is a thermodynamic ally 
unstable region in this model, therefore the curves in Fig. |7]are 
not continued to the left beyond this point. For the magnetized 
iron model, the melting curve does not cross the surface, because 
Tm > T^crit- In this case, the open circles mark the density that 
the condensed surface would have at much smaller temperatures 
T ~ 10^ K <K Tcrit. The parts of the curves to the left of the open 
circles cannot be reached in a stationary stellar envelope. 

7.4. Thermal structure of a magnetar envelope 

The results presented above have a direct application to the cal- 
culations of the thermal and mechanical structure of neutron-star 
envelopes with strong magnetic fields. Figure |8] illustrates the 
structure of a typical magnetar envelope with the ground-state 
nuclear composition. For illustration we have assumed that the 
magnetar has mass 1 .4 Mq and radius 12 km, and the considered 
patch of the stellar surface has effective temperature 10^^ K and 
magnetic field B = lO'^ G inclined at 45°. The top panel shows 
the thermal structure of the envelope, which has been calculated 
by numerical solution of the system of heat balance equations, 
taking the General Relativity e ffects and neutrino emission into 
account (Potekhin et al.ll2007l). The middle pan el presents the 
ion charge Z as function of p (iRuster et al.ll2006 V In the bottom 
panel (analogous to Fig. |6]l we plot several reduced thermody- 
namic functions of p and T along the thermal profile (i.e., taking 
T from the top panel), starting at the condensed solid surface. 
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Fig. 8. Structure of a magnetar envelope having the ground-state 
nuclear composition, the effective temperature 10^'^ K, and mag- 
netic field B = 10'^ G, inclined at 45° to the surface. Top panel: 
thermal profiles calculated using the present EOS (red solid line) 
and the EOSs where the Coulomb nonideality is either neglected 
(green dotted line) or treated without account of magnetic quan- 
tization (blue dashed line, which is superposed on the red solid 
line). The melting temperature is drawn by the oblique dot- 
dashed line. Thin vertical dot-dashed lines mark the points of 
phase transitions from solid to liquid and back to solid state. 
Asterisks mark the ends of the convective segme nt, which is 
pointed out by the arrow. Middle panel: ion charge ( Riister et al.l 
2006). Bottom panel: reduced pressure, heat capacity, and loga- 
rithmic derivatives of pressure. 



The temperature quickly grows at the solid surface and 
reaches the melting point at the depth z ~ 1 cm. Thus, at the 
given conditions, the liquid ocean of a magnetar turns out to be 
covered by a thin layer of "ice" (solid substance). We treat the 
solid crust as immobile, but the liquid layer below the "ice" is 
convective up to the depth z ~ 1 m. We treat the convective 
heat transport through this layer in the adiabatic approximation 
(ISchwarzschildlll958l) . The change of the heat-transpoit mech- 
anism from conduction to convection causes the break of the 
temperature profile at the melting point. We underline that such 
a treatment is only an approximation. In reality, the superadia- 
batic growth of temperature can lead to a hydrostatic instability 
of the shell of "ice" and eventually to its cracking and fragmenta- 
tion into turning-up "floes". This can result in transient enhance- 
ments of the thermal luminosity of magnetars. 

The temperature profile flattens with density increase, and 
the Coulomb plasma freezes again at the interface between the 
layers of *Ni and ^''Kr at p = 1.5 x lO** g cm"^ (z = 73.8 m). 
These phase transitions do not cause any substantial breaks in 
Xp, XT, or Cv/Niks, because the Coulomb plasmas have similar 
structure facto rs in the hquid and crystalline phases in the melt- 
ing region (cf. lBaiko et"al]ll998l) . 
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At the boundaries between layers composed of different 
chemical elements, the reduced thermodynamic functions do 
not exhibit substantial discontinuities, except for the abrupt in- 
creases of P/riikBT at the interfaces ''*'Ni/**''Kr (p = 1.5 x 10'' 
gcm"3) and ^^Ni/'^'^Mo (p = 1.32 x 10" g cm"^), which are 
caused by the decreases of with the large jumps of A (of 
course, the non-normalized pressure P is continuous). The spe- 
cific heat per ion Cv/Niks is almost continuous at these inter- 
faces, which means that heat capacity of unit volume abruptly 
decreases. The drop of xt at the ^**Ni/'-**Mo interface is due to 
the same decrease of ni, which leads to the decrease of the ionic 
contribution that mostly determines dP/dT at the strong degen- 
eracy. 

The oscillations of the reduced thermodynamic functions 
(most noticeable for Xp ^^ndxr) correspond to consecutive pop- 
ulation of excited Landau levels by degenerate electrons with 
density increase, analogous to the oscillations in Fig.|6] 

The magnetic effects on the nonideal part of the plasma ther- 
modynamic functions have almost no influence on the temper- 
ature profile in the magnetar envelope, as illustrated in the up- 
per panel of Fig. [8] where the corresponding solid and dashed 
lines virtually coincide. For comparison, the dotted line in the 
upper panel shows the result of a calculation totally neglecting 
the Coulomb nonideality. In this case, the profile is quite differ- 
ent at low densities, where there is no solid surface anymore. 
However, even in this case the thermal profile is almost the same 
at large p. This means that the Coulomb nonideality has a minor 
impact on the relation between the internal and effect ive temper- 
atures and therefore on the cooling curves (lYakovlev & PethicS 
l200l . but it c an be important for th e shape of the thermal spec- 
trum (cf., e.g.. |Potekhin et al.ll2012l) . 



8. Conclusions 

We have systematically reviewed analytical approximations for 
the EOS of fully ionized electron-ion plasmas in magnetic fields 
and described several improvements to the previously published 
approximations, taking nonideality due to ion-ion, electron- 
electron, and electron-ion interactions into account. The pre- 
sented formulae are applicable in a wide range of plasma pa- 
rameters, including the domains of nondegenerate and degener- 
ate, nonrelativistic and relativistic electrons, weakly and strongly 
coupled Coulomb liquids, classical and quantum Coulomb crys- 
tals. As an application, we have calculated and discussed the be- 
havior of thermodynamic functions, melting, and latent heat at 
crystallization of strongly coupled Coulomb plasmas with the 
parameters appropriate for cooling white dwarfs and envelopes 
of nonmagnetized and strongly magnetized neutron stars. We 
have also shown that a typical outer envelope of a magnetar can 
have a liquid layer beneath the solid surface. 

We have made the Fortran code that realizes the analyti- 
cal approximations described in this paper freely available at 
http : //www . iof f e . ru/astro/EIP/. 
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Appendix A: Nonideal part of the free energy of 
electrons 



lTanakaetal.l (fl98?) calculated the interaction energy of the elec- 
tron fluid at finite T and presented a fitting formula that repro- 
duced their numerical results as well as the results of other au- 
thors in various limits. Subsequ ently the behavior of the fit at 
r <K Tp was improved by Ichimaru et al.l (Il987h . The result 
reads 



2A 



2(dB + CA) ^ 



fD 



'2f ^/^^+d\ Id 

arctani - arctan — 

D D 



-(l-^jM/Fe + fl'Vr^+i), 



(A.l) 



where A^b-gd, B = a-g,C = 2- d^/f, D = ^J4f - d^, 
and a, b, d, f, and g are the following functions of = T/Tp (at 
B = 0): 

-^1 tanh - 
4n^j 

0.75 + 3.0436302 - 0.0922761^ + 1.70350^ 

^ 1 +8.31051(92+5.110504 ' 

0.341308 + 12.0708 6|2+ 1.148889 0* r , 1 

b = T — ye tanh — , 

1 + 10.495346 612-H 1.326623 04 ' 

, 0.614925+ 16.99605502+ 1.48905604 ^ 1 

d = r V0 tanh — , 

1 + 10.10935 02 + 1.2218404 ' 



/ 



0.539409 + 2.522206 02 + 0.17848404 



tanh 



1 +2.55550102 + 0.146319 04 
g = 0.872496 + 0.025248 exp(- 1/0). 

The accuracy of Eq. (lA.ll i is 1%. 

In a quantizing magnetic field, we replace the argument 
in these expressions by the quantity 0* defined by Eq. ( l74l l. as 
explained in Sect. 15.11 

Appendix B: Nonideal part of the free energy of 
ions in the rigid background 

B. 1 . Coulomb liquid 

For the reduced free energy fa = Fa jNik-^T of the classical OCP, 
we have the following analytical formula (Paper I): 



I 



.(0) 



^|f{A7+Y) - A2\n[^|f|A'2 + VTTTtaI) 



+2A3 [Vf - 



arctan 



Vf] 



+Bi 



r-B2ln|l + — 



B3 / r2 

+ — In 1 + — 
2 I B4 



(B.l) 



where Ai = -0.907347, A2 = 0.62849, B, = 0.0045, B2 = 170, 
B3 = -8.4 X 10-^ B4 = 0.0037, and A3 = - V3/2 - Axl ^fM. 
The derivative 



(0) 



d\nT 
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VfTa; 



r + 1 



Bir2 



B3r2 



Y + B2 r2 + B4 ' 



(B.2) 



reproduces Monte Carlo calculat ions of the red uced internal en- 
ergy UiilNik^T at 1 < r < 190 (ICaillolll 19991) within the accu- 
racy of these calculations, < 10"^. For any values of the coupling 
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parameter in a liquid OCP, < F < 200, the fractional eiTor of 
the approximation (IB. Il l does not exceed 2 x 10 

The classical treatment of ion motion is justified at T » Tp. 
One can extend the applicability range of the analytical EOS to 
r ~ Tp by the Wigner-Kirkwo od quantum corrections ( Wignen 
[T93I iLandau & LifshitdfToSOl) . The lowest-order coiTection to 
the reduced free energy is 

/f = 772/24. (B.3) 

The next-order coiTection oc Ji^ was obtained by 
iHansen & Vieillefo sse ( 1975|). It can be written as 



<4) 



- 2.085 X 10"^ - 



2.411 X 10- 



0.001288 



ri/2 r 

1.353 x 10-"* 0.002476 0.00276\4 



p3/2 



r2 



(B.4) 



This expression, unlike Eq. (IB.3l l, is not exact. Both corrections 
have limited applicability, because as soon as 77 becomes large, 
the Wigner expansion diverges and the plasma forms a quantum 
liquid, whose free energy is not known in an analytical form. 
Therefore we use only the lowest-order correction ([B.3) . i.e., 
/ii ~ /j^"* + /P'— In a magnetic field, Eq. ( IB.3b is replaced by 
Eq. (123 (Sect'l5:2l). 

B.2. Coulomb crystal 

The reduced free energy of an OCP in the crystalline phase is 
given by Eq. (l32l l. wh ere the first thr ee terms describe the har- 
monic lattice model (B aiko et"ar '2001). For the bcc crystal, we 
have Co = -0.895 929 255 68 and m = 0.51 1 3875, and for f± 
the following fitting formula can be used: 



/th = Xln(l 



k=l 



Airf) 



(B.5) 



where ai = 0.932446, = 0.334547, 03 = 0.265764, 
A(7^) 



1 -H 0. 1 839 77 + 0.593 586 77^ -h 0.005 48 14 77^ 



+5.01S 13 X 10-"* 77"* + 3.9247 x 10"^ 77'' 

-H5.8356 X 10"" 77^ 

B(77) = 261.66 H- 7.079 97 77^-^0.0409 484 77"* 

-H3.973 55 X 10""* 77^ -h 5.1 1 1 48 x 10"^ 77^ 

-H2.19749X 10"^ 77' + 1.866 985 x 10"''77'' 

-H2.787 72X 10"'^77". 

The Taylor expansion of Eq. (IB. 5b at small 77 is consistent with 
the Wigner correction (IB. 3b . However, the next Taylor term ~ 77^ 
is absent in the Wigner expansion, and therefore Eq. (IB. 5b does 
not reproduce higher-order Wigner corrections. Nevertheless, 
approximation ( IB. 5b is very accurate: it reproduces the numeri- 
cal results in Baiko et al. (2001) with fractional deviations within 
5 X 10 '', and its first and second derivatives reproduce the calcu- 
lated contributions to the internal energy and heat capacity with 
deviations up to several parts in 10^'. Other types of simple lat- 
tices are described by the same expressions with slightly differ- 
ent parameters (see Baiko et al. 2001). 

Anharmonic corrections for Coulomb lattices were studied 
in a number of works (see Papers 1 and 11 for references). In the 
classical regime 77 0, we have chosen one of 1 1 parametriza- 
tions proposed bv lFarouki & Hamaguchil(ll993l) : 



Ok 



(B.6) 



where fli = 10.9, a2 = 247, and 03 = 1.765x10 . A continuation 
to arbitrary 77, which is consistent with available analytical and 
numerical results for quantum crystals, reads (Paper II) 



/ah = /r(r)e-°'"''"^-o.i2 77Vr. 



(B.7) 



Superstrong magnetic fields can significantly change these 
expressions under the conditions > 1 and 77 > 1. Analytical ap- 
proximations for the free energy of a harmonic Coulomb crystal 
in quantizing magnetic fields are derived in Sect. |6] Analogous 
results for the anharmonic corrections are currently unavailable. 

Appendix C: Electron polarization corrections 

C.1. Coulomb liquid 

The screening contribution to the reduced free energy of the 
Coulomb liquid at < F < 300 has been calculated by the HNC 
technique and fitted by the expression (Paper I) 



-F, 



CdH Vrr + CTFflF^gl/ll 



(C.l) 



where cdh = (Z/ y/3)[(l + Zf'^ - I - Z^'^] ensures exact 
transition to the Debye-Hiickel limit at F — > 0, ctf = 
(18/175) (12/7r)2/3z^/3 (j _ 2-1/3 +o.2Z-i/2^ f^^^ jj^g n umerical 
data a t large F and reproduces the Thomas-Fermi limit (ISalpeterl 
1961) at Z ^ 00, the parameters a = 1.11 Z" "*", b ^ 0.2 + 
0.078 (lnZ)2, and v = 1.16 H- 0.08 InZ provide a low-order ap- 
proximation to Fig for intermediate and F. Finally, the func- 
tions 

gi ^ 1+ 0.78 [21 + Fe(Z/r,)3]"' (rjzy'\ 

^2 = 1 + — 11 + 



0.001 Z2H-2Fe/ 1 -i-6rr 
improve the fit at relatively large r,, and 

^ ^ 1 + 

' 1 + 0.18Z-i/4xr + 0.37Z-i/2x? + x?/5 

is the relativistic correction, as well as y^' in the denominator. 
C.2. Coulomb crystal 

The screening contribution to the reduced free energy of the 
Coulomb crystals was evaluated using the semiclassical pertur- 
bation approach with an effective structure factor (Paper I) and 
fitted by the expression (Paper II) 



/ie = -fMr[i +A(^,)[e(77)/r]'}, 

where 



(C.2) 



/oo(jc) = ajpZ^'^bi y/l+b2/x^, 

^ ^73 + 17.9x2 

A(x) = — ——, 

1 + b4X^ 

2(77) = [in (1 + [in _ _ 2)e-^'-''"''Y"^ ' 

the parameter otf = 0.00352 is related to cjp in Eq. ( IClb . and 
parameters s and b\-bn depend on Z: 

s = [iH- 0.01 (lnZ)^/2-H 0.097 Z"2]"', 
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= 1 - 1.1 866 Z"°-2''^ + 0.27 Z"', 
, 2.25 1 + 0.684 +0.222 Z'' 

09 = 1 H 7T7 ; , 

Zi/3 1 +0.222 Z6 

/?3 = 41.5/(1 +lnZ), 
bA = 0.395 lnZ + 0.347 Z^^^l 

Here, the numerical parameters are given for the bcc crystal; 
their values for the fee lattice are slightly different (Paper I). 
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